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The response of glasses to mechanical loading often leads to the formation of 
inhomogeneous flow patterns. Among them, shear bands, associated with strain 
localization in form of band-like structures, are ubiquitous in a wide variety of mate¬ 
rials, ranging from soft matter systems to metallic alloys. They can be precursors to 
catastrophic failure, implying that a better understanding of the underlying mecha¬ 
nisms of shear banding could lead to the design of smarter materials. Here, molecular 
dynamics simulations are used to reveal the formation of shear bands in a binary 
Lennard-Jones glass, subject to a constant strain rate. At a critical strain, this sys¬ 
tem exhibits for all considered strain rates a transition towards the formation of a 
percolating cluster of mobile regions. We give evidence that this transition belongs 
to the universality class of directed percolation. Only at low shear rates, the per¬ 
colating cluster evolves into a transient (but long-lived) shear band with a diffusive 
growth of its width. 

Introduction. An external shear field leads in general to a rejuvenation of the glass state, 
transforming the amorphous solid into a flowing fluid una. Under a constant strain rate, 
the transition to plastic flow can be located via the dependence of the shear stress on the 
applied strain. It is marked by a maximum in the stress-strain relation which for a simple 
planar Couette flow geometry occurs typically at a strain of the order of 0.1 [2HS]. Beyond 
this maximum, the system evolves into a steady-state regime where the system displays a 
homogeneous flow pattern, e.g. in the case of planar Couette flow manifested as a linear 
velocity profile. However, the strain necessary to reach this steady-state regime depends 
on many factors such as the history of the initial undeformed glass state and the applied 
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strain rate mm- In the transient regime before the steady state is reached, the occurrence of 
spatially inhomogeneous flow patterns is very common. Often, such inhomogeneous response 
leads to the formation of shear bands, with co-existing regions of contrasting mobilities 
spanning large scales [8]. Even when the mechanical loading is switched off, the shear-band 
regions can be imprinted as frozen-in structures in the resulting solid and thus strongly affect 
its material properties [91 E]. Therefore, the microscopic understanding of the formation 
and growth of shear bands in glassy solids is a highly-debated issue in material science. 

Despite the fact that shear banding is ubiquitous in amorphous materials [HlinHIl], the 
microscopic processes that lead to the formation of such complex structures still remain ill 
understood. Unlike the case of micellar systems or granular materials, where shear bands 
mainly form in the regime of large applied shear rates, driven by a coupling between struc¬ 
tures and external shear [HHIZ], the scenario in glasses seems to be different. The inhomo¬ 
geneous response in the case of amorphous solids tends to be observed at small shear rates, 
with the spatio-temporal extent becoming more prominent as one approaches the yielding 
threshold EEllIHl fT9] . And, so far, there is no indication of any underlying macroscopic 
mechanical instability driving the formation of shear bands. While there have been several 
observations in experiments [2011^ . numerical simulations [221126] and phenomenological 
models EZIEH], a quantitative analysis that allows to predict the conditions under which 
shear bands form in glasses is lacking. 

The question about the origin of shear bands is intimately related to the question how 
flow is initiated in an amorphous solid under applied shear. Recent studies suggest that the 
response of the glass to the applied strain is governed by local heterogeneities that are either 
already present in the undeformed solid or form during the initial application of the shear 
held [H Elll] . These heterogenities are associated with “hot spots” of higher mobility that 
grow while the strain of the system increases m [30] . For athermal amorphous systems, 
evidence has been given that the transition towards plastic how is provided by an avalanching 
of the mobile hot spots [OTHOO] . Imprints of that have also been reported in thermal systems 
[on 100] . In fact, such avalanching is well-known in the context of self-organized criticality 
[06] . A prominent example is the sandpile model [36]. Recently, it has been shown that a 
driven version of this model belongs to the directed percolation (DP) universality class [07] . 
In fact, as put forward by the DP conjecture [OH [39], DP universality is believed to be very 
robust, applying to a broad class of non-equilibrium phase transitions [iQ]. 
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Thus - although not addressed so far - it is a natural question whether the onset of flow 
in a sheared glass is linked to a DP transition. In this work, we show that this indeed the 
case and thereby we elucidate the microscopic mechanism that leads to the formation and 
growth of shear bands in glassy systems. 

Using a model glass forming system, we use large-scale molecular dynamics simulations 
to study the mechanical response of a quiescent amorphous solid when an external shear rate 
is imposed on it. We track the locations of regions of large mobilities and thereby we reveal 
the existence of the DP transition driven in the direction of applied shear. While the DP 
transition is seen for all shear-rates, visible shearbanding only only for small shear rates near 
the yielding threshold. We demonstrate that the glass eventually fluidizes by the long-time 
diffusive invasion of the shear band into the rest of the system, with the diffusion timescales 
dependent on the imposed shear rates. Thus, we provide a quantitative description of how 
initial local mobilities build up to the eventual macroscopic flow of the material. 

Results. When the shear is applied to the quiescent glass at t = 0, the material deforms 
exhibiting the typical stress (a) vs. strain (jt) response [shown in Fig. [^a) for different 
imposed shear rates 7 ] with the height of the overshoot depending on 7 IS]- The measured 
steady-state stress as a function of the imposed shear rate is shown in the inset of Fig. [^a); 
it has the typical Herschel-Bulkley form [B]. The corresponding single particle dynamics, 
during this onset of flow, can be quantified by measuring the non-affine mean-squared dis¬ 
placement (MSD), Ar^^, in the direction transverse to the applied shear; the data is shown 
in Fig. [^b). The particles undergo ballistic motion at early times and are then caged, before 
the occurrence of a super-diffusive regime, prior to diffusion. The onset of super-diffusion 
occurs around the stress overshoot in the stress-strain curve, when the built-up stress is 
released via the particles breaking their local cages to subsequently diffuse [5|. 

In experiments, flow heterogeneities are often diagnosed via the spatial prohles of lo¬ 
cal velocities imiig. Similarly, we measure the spatial profiles of the local flow velocities 
(averaged over strain intervals of 0.5%), V 3 ;{z), for an imposed shear rate of 10“^, at differ¬ 
ent times after the imposition of shear (marked on the corresponding stress-strain curve in 
Fig-i a)). The spatial profiles are shown in Fig. [^c) for one of the initial states in our en¬ 
semble. In the elastic regime, at jt = 0.02, the velocity profile is linear, but starts deviating 
from this shape as the stress overshoot is approached at yt = 0.06. This deviation becomes 
stronger in the transient regime (at jt = 0.5, shown in blue diamonds). Interestingly, the 
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FIG. 1. (a) Stress-strain response of the glass at temperature T = 0.2 for an age of = 10“^ 
and sheared with a constant strain rate. Data shown for different shear-rates: 7 = 10“^, 3 x 
10“^ and 10“^. The inset shows the flow curve. Red solid line shows the fitting with Herschel- 
Bulkley form — 0-3974 -|- 2.29637*^'^^. (b) Variation of z-component of the MSD of large 

particles with strain. Brown dotted line marks ^th{= 0.02). (c) Velocity profiles, for 7 = 10“^, at 
five different strain values marked in (a). Green solid line represents the expected linear profile, 
(d) MSD maps at 'jt = 0.06,0.5 for 7 = 10“^. (e) Maps of local strain corresponding to MSD maps 
shown in (d). (f) Maps of local mobility corresponding to (d). Mobile regions are marked in blue 
while immobile regions are marked in green. 

velocity profile regains its linear shape as plastic flow sets in (at jt = 1 . 0 , shown in orange 
stars). The observation of increased heterogeneity, after the stress-overshoot, is consistent 
with experimental observations [ 12 ]. 

However, local velocity prohles only captnre the short-time heterogeneities in dynamics. 
In order to obtain a more cnmnlative pictnre from t = 0, we look at the spatially resolved 
maps of Arz‘^{t) [23] HI], by dividing the system into small cells (see Methods for more 
details). In Fig. [^d), we show the time-evolntion of snch a map, for an initial state nnder 
the imposed shear rate of 7 = 10 “^. At a strain of jt = 0.06, the local dynamics is nearly 
homogeneons on this scale. However, at jt = 0.5, spatially heterogeneons dynamics is 
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FIG. 2. Occurrence of percolation transition, (a) Variation of Pspan with p, in the box of dimension 
20 X 20 X 80, for strain rates 7 = 10“^, 10“^, 3 x 10“^, 10“^, 3 x 10“® and 10“^. (b) Percolating 
cluster at the critical point for the different strain rates: 7 = 10 “^, 10 “^ and 10 “^ (left to right), 
(c) Variation of p with strain for all strain rates shown in (a). Orange dotted line corresponds to 
the 'jt = 0.5. (d) 2-Dimensional slice of MSD map for trajectories shown in (b), at jt = 0.5. 


observed, with the more mobile particles localised in a shear-band-like structure spanning 
the X — y plane. We can also construct similar maps of local strain [szx = dAvz/dx), which 
exhibit a localisation behaviour similar to the MSD maps; see Fig. [^e). Thus, large local 
MSDs are also regions of large strains. Henceforth, we use to analyse local dynamical 
properties. 

In order to quantify and characterise the spatio-temporal evolution of the mobile regions, 
we dehne a region to be mobile or not, by setting a threshold pth = 0.02 on the local Ar^^. 
As marked by the dashed line in Fig. [^c), such a choice of pth is larger than the plateau 
value in the MSD and thus corresponds to motions beyond cage-breaking [50]. We then 
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define the local mobility 'ifj as 



0 otherwise, 


1 if /i ^ 


( 1 ) 


where /i is the average MSD of particles in a sub-box. Following this convention, we digitize 
the whole system into mobile and immobile regions. The mobility maps corresponding to 


Fig. [^d) are shown in Fig. [^f). 


We now demonstrate that a percolation transition occurs with increasing strain, involving 
these mobile regions. As the system evolves under the applied shear rate, we monitor the 
fraction of mobile cells, p, at any given instant. Figure[^d) suggests that such mobile regions 
do form clusters. Thus, we compute what fraction of these mobile regions, Pspan, is part of 
a cluster that spans the system, recalling that such an observable is the order parameter 
for determining the occurrence of a percolation transition. In Fig. [^a), we plot Pspan as a 
function of p, which shows that beyond a critical fraction p^, all the mobile cells are part of 
such a spanning cluster. This indicates the occurrence of a percolation transition of these 
mobile cells. Furthermore, we observe that the variation of Pc is nearly independent of the 
imposed shear-rate, as seen in Fig. |^a) for a wide range of 7 . Thus, the percolation process 
is generic to the system’s response under shear. In Fig. [^b), we visualise the spanning 
clusters for three different shear rates, at the Pc corresponding to each 7 , starting from the 
same initial state. Here, we note that the spanning across the system-size occurs in the 
direction of flow, which we discuss further later. 

Next, for different imposed 7 , if we monitor the numerical growth of mobile cells with 
increasing strain, a distinct variation is revealed; see Fig. |^c). For example, at a strain of 
■yt = 0.5 (marked by orange dotted line), we see that around 50% of the sites are mobile 
at low strain rates while the same is close to 75% at high strain rates. This implies that 
subsequent to the percolation transition (at the critical strain corresponding to pc, which 
is around 0.0685 for 7 = 10“^), the spatial heterogeneity of activity is more long-lived for 
smaller strain rates, as is shown in Fig. |^c), via cuts in the x — z plane of the local MSD 
maps corresponding to the evolving trajectories of the states shown in Fig. |^b). While for 
the largest shear-rate, mobile regions proliferate in the system, for the smaller shear-rate it 
is more localised and takes the form of a well-structured shear band. Thus, one can infer 
that the local dynamics, post-percolation, changes with decreasing shear-rate. 
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FIG. 3. Occurrence of directed percolation (DP) transition, (a) Variation of Pspan with p for four 
different system sizes 10 x 10 x 40 (black circles), 20 x 20 x 80 (magenta squares), 40 x 40 x 40 (green 
diamonds), 50 x 50 x 50 (blue triangles) for 7 = 10“^. (b) Finite size scaling of critical points for 
the systems sizes shown in (a), the red solid line shows the fit according to the scaling law Pc{L) = 
0.30339539 — 2.412L“^/^'^°® corresponding to DP. Green dash-dotted line shows the same scaling 
function with parameters corresponding to standard percolation: critical point Pc{oo) = 3.116 and 
V = 0.8765. The blue dotted line shows the asymptotic Pc{oo) for DP which is 0.30339539. (c) 
Gluster size distribution around pc for the system sizes shown in (a), data shows a power-law decay 
with exponent —2.4 which is consistent with the prediction for DP. 


To clarify the nature of the percolation transition, we determine the critical point for 
the percolation process, Pc, using finite size scaling. In Fig. [^a), we show how Pspan varies 
with p for four different system sizes - we observe that the onset of percolation shifts to 
larger values of p with increasing system size. By fitting the obtained threshold for different 
system sizes, using the finite size scaling function Pc{L) = Pc{oo) + (see Fig. |^b), 

red line), we obtain an estimate of Pc(oo) = 0.30339539 and the exponent i/y = 1.106, which 
corresponds to a DP transition |13]. To compare, the corresponding numbers for standard 
percolation are Pc{oo) (= 0.3116) and z/ (= 0.8765); the finite-size scaling function using 
these parameters does not fit our data (see Fig. [^b), green line). Furthermore, we compute 
the size distribution of clusters of mobile cells in the vicinity of the percolation transition 
(Fig. i c)). As expected, the distribution has a power-law shape, with the corresponding 
DP exponent of 2.39 well characterising the distribution gaisi. Thus, the percolation of 
the active regions, in this regime of flow, is a directed process, driven in the direction of the 
external shear. 

Thus far, we have discussed the existence of an underlying percolation process of the 







FIG. 4. (a) Two dimensional projections of MSD maps showing the evolution of shear band with 

time, for imposed 7 = 10“'^, at strains of jt = 0.2, 0.5,1.0, 3.0. (b) Variation of width of the band 
^b, scaled by the box length Lz, with 'jt = 10“^, 10“^, 3 x 10“^ and 10“®. The solid red line 
corresponds to 


mobile regions, which occurs for the entire range of 7 that we have studied. Now, we will 
focus on how the dynamics proceeds once the percolating cluster has formed. In order to 
quantify that, we construct two-dimensional projections of the MSD maps. In Fig. |^a), 
we show the time evolution of the local mobility, via such maps, for 7 = 10“^. We observe 
that the 2 ;-width of the shear band increases with time and eventually the entire system 
becomes mobilised. By locating contiguous layers of mobility, we identify the shear band 
and, thereafter, by marking the interfaces of this band, we measure how the band-width, 
.^b, evolves with time. For different imposed shear rates, this time evolution is shown in 
Fig. I^b). We see that initially grows quickly and then eventually it reaches a regime 
where the data can be htted with ,^b ~ implying that the propagating interface of the 
shear band has a diffusive motion. The diffusion constant is dependent on the imposed 
shear rate. The smaller the shear rate 7 , the slower is the diffusion, which leads to more 
long-lived heterogeneities, as discussed earlier. For the largest shear rate shown, 7 = 10“^, 
the diffusive regime is very short-lived as the band quickly spans the entire system. For 
even larger shear rates, beyond percolation, mobile regions quickly appear everywhere and 
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fluidizes the whole system, and, as a consequence, a shear band is not clearly discernible. 

Conclusions. In this work, we have explored how a model glass, subject to a constant 
strain rate, evolves from a quiescent state to plastic flow. We have shown that this process is 
initiated by a DP transition. Under shear, hot spots form in the amorphous solid, i.e. local 
regions in the system where particles have undergone large non-affine displacements. These 
are local structural changes transforming the initial quiescent glassy state which thus cannot 
be regained via thermal fluctuations. Unlike liquids, this is an essential feature of the 
glassy state, where a strong non-linear response to the external shear over-rides thermal 
fluctuations. Under continuing deformation, such hot spots percolate in the direction of the 
applied drive. Such a scenario conforms to the DP conjecture [3H1 ES], whereby a system, 
with no quenched disorders, transforms from a fluctuating state into an absorbing state. 
In our case, the dominance of the external shear leads to the system getting irreversibly 
trapped in an absorbed state of the percolating cluster of hot spots. 

The subsequent growth of the active regions after the DP transition depends strongly on 
the applied shear rate. At larger shear rates, we observe a quick proliferation of mobile spots 
leading to a fluidisation. On the other hand, the process is asymptotically much slower at 
smaller shear-rates, with the mobile front slowly diffusing into the rest of the material. Thus, 
the existence of an initial fast timescale for the burst in mobility in the flow direction along 
the neutral plane and the subsequent slow timescale for the spread in the shear-gradient 
direction, leads to the sustenance of shear bands with long timescales for small shear rates. 


I. METHODS 


We consider a binary mixture of Lennard-Jones (LJ) particles (say A and B) with 80:20 
ratio. This is a well-studied glass former. Particles interact via LJ potential which is dehned 
as: 


= (l>ap{r) - (j)afi{Rc) “ - Rc) 

= 4e„/3 [{crafslrf^ - {cro^p/rf] r < R 


r=Rc 

Cl 


( 2 ) 


where a,/3 = A, B. Interaction among particles is dehned as caa = 1-0, cab = I-Scaa, 
cbb = O.Scaa- Range of interactions is given as ctaa = 1-0, ctab = O.Sctaa, o'bb = 0.88ctaa 
and Rc = 2.5(Jaa- Masses of both type of particles are equal, i.e., tua = All 
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quantities are expressed in LJ units in which the unit of length is ctaa, energy is expressed 
in the units of caa and the unit of time is -^/mcr^^/eAA- More details about the model and 
parameters can be found in Ref. |15] . 

We perform molecular dynamics (MD) simulation in the NVT ensemble using the package 
LAMMPS (“Large-scale Atomic/Molecular Massively Parallel Simulator”) [17j. Different 
geometries are considered, placing the particles in boxes of dimensions 10 x 10 x 40, 20 x 
20 X 80, 40 X 40 X 40 and 50 x 50 x 50. Temperature is maintained via dissipative particle 
dynamics (DPD) thermostat [48] . 

Our method for the preparation of glass is as follows: First we equilibrate high temper¬ 
ature initial conhguration at a temperature T = 0.45 in the super-cooled regime and then 
quench it to a temperature T = 0.2 below the mode coupling transition temperature [46] . 
We wait until = lO'^ and apply shear on x-z plane in the direction of x with different 
constant strain rates 7 = 10 “^, 10 ’ -^3xl0-^10-^ , 3 X 10 ® and 10 To simulate a sheared 
bulk glass, we use Lees-Edwards periodic boundary conditions [4U] . 

Maps of local MSD. To get the MSD maps, we divide the simulation box into small 
cubic sub-boxes having linear size of ctaa- At any time t, we calculate the average MSD of 
the particles populating each sub-box at f = 0 (unsheared glassy state). As discussed earlier, 
we measure the ^-component of MSD of each particle. The corresponding three-dimensional 
map, thus constructed, is shown in the different plots. 

Maps of local strain. To plot the strain maps, again we divide the simulation box into 
small cubic sub-boxes, as in the case of MSD maps. The numerical derivative of x-component 
of displacement with respect to x, i.e. e^x = dAvz/dx is calculated. This derivative, which 
gives strain Czx, is plotted for each sub-box to construct the map. 

Identifying the interface of shear-band. To characterize shear bands, we choose a 
slightly higher threshold on MSD. In the present work, it is taken as pth = 0.1. We divide 
the simulation box into x-y layers of thickness one particle diameter and calculate y, which 
is the x-component of MSD, for each layer. We assign each layer a value V’ =1/0 depending 
on whether /i is greater/lower than yth- To get the size of shear band we count the number 
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of adjoining layers for which '0 = 1. 
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